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Abstract: The movement of many organisms can be described as a random walk at either or both the individual and 
population level. The rules for this random walk are based on complex biological processes and it may be difficult to 
develop a tractable, quantitatively-accurate, individual-level model. However, important problems in areas ranging from 
ecology to medicine involve large collections of individuals, and a further intellectual challenge is to model population- 
level behavior based on a detailed individual-level model. Because of the large number of interacting individuals and 
because the individual-level model is complex, classical direct Monte Carlo simulations can be very slow, and often of little 
practical use. In this case, an equation-free approach |24| may provide effective methods for the analysis and simulation of 
individual-based models. In this paper we analyze equation-free coarse projective integration. For analytical purposes, we 
start with known partial differential equations describing biological random walks and we study the projective integration 
of these equations. In particular, we illustrate how to accelerate explicit numerical methods for solving these equations. 
Then we present illustrative kinetic Monte Carlo simulations of these random walks and show a decrease in computational 
time by as much as a factor of a thousand can be obtained by exploiting the ideas developed by analysis of the closed 
form PDEs. The illustrative biological example here is chemotaxis, but it could be any random walker which biases its 
movement in response to environmental cues. 



1 Introduction 

In current complex systems modeling practice, we are often presented with a model at a fine level of description 
(atomistic, stochastic, individual-based), while we want to study the behavior at a macroscopic coarse-grained 
(continuum, population) level. This situation frequently arises in the modeling of biological dispersal, where 
significant progress is being made in modeling at the individual organism/cell level, while the derivation of 
the corresponding closed, macroscopic population level equations remains very difficult, and lags far behind in 
development. The example here is bacterial chemotaxis, for which much is known about signal transduction and 
motor behavior of individual cells, but only in a limited number of cases can one rigorously derive equations 
describing behavior of bacterial populations |1U1 ITT] . Usually one can develop a suitable cell-based stochastic 
model, and would like to obtain population-level information without having a coarse-grained evolution equation. 
Computational methods for obtaining an approximation to the macroscopic evolution without explicitly obtaining 
equations have been developed |241 1151 ITY1 137| . The main idea is to use short bursts of appropriately- initialized 
computations using the detailed, fine-scale model, followed by processing of the results to obtain estimates of 
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the desired macroscopic quantities such as the spatial distribution of the number density, time derivatives, and 
various measures of the sensitivity of the solution with respect to parameters. 

The first, and probably most important step of this equation-free approach is to determine what are the 
appropriate variables in terms of which one could hope to close macroscopic evolution equations. Typically 
these variables are a few slowly-evolving lower moments of the many-particle distribution function (e.g., cell 
density for chemotactic movement |10| , species concentrations for reaction-diffusion problems |14j , or density and 
momentum fields, the zeroth and first moments of the distribution of molecules in velocity space, for the Navier 
Stokes equations HOE!)- In most cases, knowledge of the level of closure (the number and identity of variables with 
which one can write a deterministic model for the process) comes from extensive experimental experience and 
observation, long before it is rigorously justified by theory. In the equation-free approach, the simplest conceptual 
path for selecting the appropriate observables as state variables is to perform a homotopy between conditions at 
which an accurate closed equation is known and validated, and the conditions of interest, when this is possible. 
Model reduction in general is based on the assumption that, after rapid initial transients, higher order moments 
of the evolving distributions can be approximately-represented as functionals of the slow, "master" ones - the ones 
in terms of which we write the closed equations. The closure is thus embodied in a "slow manifold": a graph of 
a function (in moment space) which, given the values of the few governing lower moments, provides the "slaved" 
higher order moment values. Separation of time scales between the rapid equilibration of the higher, slaved 
moments, and the slow evolution of the "master" ones underpins the derivation of closed, reduced, macroscopic 
equations. The idea then is to design computational experiments with the fine scale simulator that test for this 
separation of time scales, and suggest variables capable of parametrizing the slow manifold |371 130j . 

As is discussed in more detail in [37], it is possible, using matrix- free iterative linear algebra methods, to 
estimate, using direct simulation, characteristic relaxation time scales for the problem (at least in the neighborhood 
of a particular equilibrium). These time scales, and the eigendirections corresponding to them, can guide the 
modeler in deciding the number, and even the selection of variables capable of parametrizing (at least locally) 
this slow manifold. In this process, homotopy and knowledge of the appropriate parametrizing variables in some 
region of operating parameter space, gives us a starting point for variable selection. In the more general, and 
much more difficult case in which we begin a completely new problem and have no initial knowledge of what 
might be good "order parameters" in terms of which to attempt to close macroscopic equations, the alternative 
is to use data processing techniques on extensive experimental (or computational experimental) runs, to try and 
develop a reasonable reduction hypothesis. Algorithms for data compression, from the more traditional principal 
component analysis to the more modern sparse kernel and diffusion map feature analysis may be useful here 
[381 I3()| . This is, however, a separate and active research subject in itself, and we will not pursue here. 

In this paper we will assume that we have enough knowledge of the problem to identify a set of variables in 
terms of which to write a closed equation. In that spirit we study the coarse integration of simple models for 
chemotaxis of cells, and we assume that the slow dynamics of the system are parametrized by cellular density. The 
main goal is to illustrate the computational gain of equation-free methods, by which we mean a large speed up of 
the stochastic simulation for a class of biologically- motivated problems involving slow dispersal of organisms/cells. 

The paper is organized as follows. In Sectional we present a brief overview of equation-free methods with 
emphasis on coarse projective integration. We present the main strategy which we will use for the analysis of coarse 
integration - namely the deterministic projective integration of partial differential equations (PDEs). Moreover, 
we show how the results of this paper can be interpreted in terms of equation-free coarse projective integration for 
kinetic Monte Carlo (kMC) simulations of random walks; and we define the gain of these methods. In Section 
we present partial differential equations modeling the dispersal of cells, and we provide two biological motivations 
of the chemotaxis system studied later. We also discuss the main mathematical properties of these equations. 
Finally, we introduce a test family of spatial signal profiles which are used in the computational examples in 
Sections 3, 4 and 5. In Section^] we study the efficiency of projective integration for different discretizations of 
the macroscopic PDE equations. We obtain a measure of efficiency (gain) of the method for different choices of 
the "inner integrator". We demonstrate a stable, signal-independent method, i.e., a method which has the same 
gain for all mathematically admissible environmental changes. We also study more accurate inner integrators, 
for which the efficiency depends on the size of the environmental signal (concentration gradients). Section R31 
contains illustrative numerical results; here we provide computations illustrating the analysis in Section^and give 
examples for which the method leads to a significant reduction in the computational time required. In Section 
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[S] we return to the original random walk problem. We discuss the application of our approach to accelerating 
the Monte Carlo simulations and present a case in which the computational time is reduced by a factor of 10 3 . 
Finally, in Section[S]we summarize the results, and mention significant generalizations. We conclude by reiterating 
the main elements of the equation-free approach as a "wrapper" around a usually slow, cell- or organism-based 
stochastic simulator, aimed at assisting in the efficient study of emergent, population-level behavior. 



2 Equation free methods - coarse integration 

Consider a large collection of randomly-walking individuals for which we have a microscopic model, and suppose 
that we want to know the time evolution of the macroscopic density N of the individuals. One approach is to 
derive partial differential equation(s) for macroscopic observables, such as the density N, and then compute the 
solution of the PDE(s) using standard numerical methods. This entails a choice of algorithm, a time step At, 
and a routine which computes the density N(t + At) from the density N(t). 

If explicit macroscopic equations are not available, we can still compute the density of individuals at time 
t + At from the density of individuals at time t using Monte Carlo simulation of the microscopic model. This can 
be done as follows. 

(a) Given the macroscopic initial density N(t), construct consistent microscopic initial conditions (initialize 
each individual so that the density is N(t)). 

(b) Evolve the system using the microscopic Monte Carlo simulator for time At. 

(c) Compute the density of individuals N(t + At) from the microscopic data at time t + At. 

Steps (a) - (c) provide an alternative path to computing N(t + At) from N(t) as illustrated in Figure^ The main 
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Figure 1: Schematic of microscopic timestepper (a) - (c). 



goal is to compute the long time evolution of cellular density N, and to that end we could simply use step (b) 
many times, i.e., we could in principle run the microscopic simulator only. However, since the biological models 
are often complex, step (b) can be very computationally intensive. Thus the key constraint is that we are in fact 
able to run the microscopic simulator only for short times. Since, we seek the long time evolution, we have to 
combine (a) - (c) with another step which can be formulated in many ways, e.g., 

(d) Using the macroscopic data for N computed in (a) - (c) , estimate the time derivative O^f {p + At) . 
Because of fluctuations due to the stochastic nature of the simulation, we may require several independent 
microscopic realizations of N(t) in part (a) to be able to accurately estimate the expected density N(t + At) 
and its time derivative. We then take advantage of the assumed smoothness (in time) of the trajectory of the 
unavailable macroscopic evolution, and take a large projective step by estimating the density N(t + At + T) 
for some T > as 

ON 

N(t + At + T) re Nit + At) + T—(t + At). 

at 

The density N(t + At + T) is then used as a new initial condition in (a). 

The algorithm (a) - (d) is called coarse projective integration, specifically, coarse projective forward Euler, and 
it can be formulated in many ways |151 1171 01) . For example we can use different methods to estimate the time 
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derivative of N in (d), or we can extrapolate using other macroscopic variables in part (d), e.g., with flux profiles 
as opposed to density profiles. In any case, the actual projective step is performed on some spatial discretization of 
these macroscopic variable profiles e.g., finite difference, finite element, or spectral decompositions of the profiles. 

The algorithm (a) - (d) can speed up the computations provided that we can safely (in terms of stability and 
accuracy) choose T ^> At and provided that the so called "lift-run-restrict procedure" (a) - (c) does not require 
excessive computation to estimate the time derivative of N. In particular, the more time is spent in part (a) - (c) 
of the algorithm, the larger T/ At in part (d) must be chosen to have the potential for computational gain. Since 
we also study modifications of (a) - (d), we will define the gain Q of the coarse projective integration method as 



Q = 



time to compute the evolution of the system by running a Monte Carlo simulator only 
time to compute the evolution of the system by coarse projective integration (a) - (d) 



(2.1) 



For example, if we need k realizations of the Monte Carlo evolution in steps (a) - (c) to compute the evolution 
of the system in the interval [t, t + At], and if we assume that the computational time of step (d) is negligible, 
then the gain Q can be simply estimated as Q ~ T ^^ t ■ On the other hand, one might argue that scaling by k 
may be too severe, since the equation we are evolving is not one for the single fluctuating realization, but for the 
expected density profile estimated, for example, as the average of k copies. 

As a first, illustrative step in the analysis of the gain of coarse integration, we will replace the stochastic part 
(a) - (c) by a deterministic operator as shown in Figure This means that we assume that we do know, at 
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Figure 2: We first analyse (i) - (ii) in parameter regimes where macroscopic equations are available. 



least for some parameter regime, a closed macroscopic equation for the expected density profile of the particular 
kinetic Monte Carlo simulation. We then replace steps (a) - (c) by a short deterministic integration (i). We run 
this deterministic integrator only for a short time At, and process its results to obtain an extrapolation in time 
(ii); we then repeat the process. In this deterministic setup, we can more easily study the dependence of the gain 
on the parameters of the model and, in particular, the gap between slow and fast eigenvalues in the spectrum of 
the equation. Assuming that most of the computational time is spent in part (i), we can rewrite the definition 
(|2.1|) in the deterministic setting (i) - (ii) as follows 



T + At 
At ■ 



(2.2) 



In the following section we introduce biologically-motivated problems for which the corresponding macroscopic 
equations are known for some, or for all, parameter regimes. 

Finally, let us mention that step (a) requires that initialization of all system variables be done consistently with 
the density profile N(t). This means that we initialize all individuals in such a way that the macroscopic density 
is equal to N(t). There are many ways to accomplish this. Ideally, we would like to initialize the remaining 
macroscopic system observables (e.g. higher moments of the cell distribution function than density, the 0-th 
moment) on a slow manifold parametrized by the density profile N(t) - that is, we would like to initialize them 
"slaved 1 to" N(t). One possible procedure (which we use in Section |5J) is schematically illustrated in Figure |31 
Here we make an initial guess for other state variables and we run the Monte Carlo simulator for a short time 

1 The underlying idea is that the set of moments of the cell distribution constitutes a singularly perturbed system, characterized 
by time scale separation: higher order moments are assumed to quickly become functionals of the low, slow, governing ones, like 
density (i.e. they quickly approach a slow manifold parametrized by density). The same key assumption also underlies the analytical 
derivation of coarse-grained, macroscopic equations. 
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Figure 3: Initialization of other state variables in step (a). 

only, then we reset the position of each individual to its initial value, keeping all other state variables unchanged. 
Repeating this procedure several times, we can find the initial condition close to the slow manifold of the system 

3 Chemotaxis 

Many organisms that move in a random walk respond to environmental signals by biasing their rules of movement. 
If we consider chemical signals in the environment, the corresponding motility behavior is called chemotaxis or 
chemokinesis, depending on whether the organism senses the direction of signal gradients directly, or responds by 
changing its speed or the frequency of turning. We will not distinguish between different terminologies, and we will 
call chemotaxis any alteration of behavior caused by the environmental cues; chemotaxis will be the illustrative 
biological example in this paper. At the population level, chemotaxis can lead to aggregation, travelling waves 
and pattern formation (see e.g., |Hj for E. coli, 0Ej for Dictyostelium discoideum) and an important task is to 
explain population- level behavior in terms of individual-based models. To do that, equation free methods may be 
suitable However our purpose here is to use the strategy for analysis outlined in Figure [SJ and to this end we 
choose a chemotactic example for which the macroscopic equations are known. First, in Section 13.11 we describe 
the simplified model of bacterial chemotaxis for which the macroscopic equations were derived in some parameter 
regimes jlOlll 1| . Next, in Section f3.2l we present an even simpler random walk, which involves directional sensing, 
and is more suitable for modeling of certain eukaryotic organisms. Here the macroscopic equations can be derived 
for any choice of parameters. These equations have the same structure as in the bacterial case. In Section 0] we 
report the results of projective integration of the chemotaxis equations. 

3.1 Bacterial chemotaxis 

Flagellated bacteria, the best studied of which is E.coli, have two different modes of motile behavior that are 
determined by the rotation of their flagella. When rotated counterclockwise, the flagella coalesce into a propulsive 
bundle that produces a relatively straight "run". When rotated clockwise they fly apart and the bacterium 
"tumbles" without significant translocation. Hence, a bacterium runs at a constant velocity for a random length 
of time, then tumbles for a random length of time, chooses a new direction at random, and repeats the process. In 
order to find food or avoid noxious substances, a bacterium increases its runs in favorable directions and decreases 
them when going in an unfavorable direction. The run length is controlled by a complex biochemical network 
[21 EH1 that involves signal transduction and alteration of an intracellular protein called CheY that controls the 
direction of rotation of the flagellar motors, and consequently changes the movement of the bacterium. 

In the absence of an extracellular signal the duration of both runs and tumbles are exponentially distributed, 
with means of 1 s and 10 _1 s, respectively and in a gradient of attractant the cell increases or decreases the 
run time according as it moves in a favorable or unfavorable direction. Since the tumbling time is small compared 
to the typical running time, we can decribe the motion of E. coli as a velocity jump process |31j . which means that 
a bacterium runs in some direction and, at random instants of time changes its velocity according to a Poisson 
process with mean turning rate 7. The turning rate is altered by CheY so we can write 7 = 7(2/1) where y\ 
denotes the concentration of the phosphorylated form of CheY. 
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Let y = (yi, y 2 , . . . , y m ) € ]R m denote the intracellular variables, which can include the concentration of 
proteins, receptors, etc., and let S(x, t) = (Si, S 2 , ■ ■ ■ , Sm) S K M denote the signals in the environment. Then 
existing deterministic models of bacterial signal transduction pathways can be cast in the form of a system of 
ordinary differential equations that describe the evolution of the intracellular state, forced by the extracellular 
signal. Thus 

|r = /(y>s) (3.1) 

where / : K m x R M — > R m describes the particular model. The equation (|3.1[1 is integrated along the trajectory 
of each cell, and the yi component of the solution together with 7 = 7(2/1) defines the random walk of each 
bacterium. 

As was noted in |l(JHllj . existing models for signal transduction and models of flagellar motor behavior involve 
tens of chemical species, which makes the problem very complicated for analysis. However, the essential aspects 
of the dynamics can be captured by a much simpler "cartoon" model which involves just two variables. For the 
"cartoon" model, one can derive closed macroscopic equations for some parameter regimes (see ^U] in ID, see 
[TT] in 2D/3D). In |1UM11| . equation 1|3.1[) and equation for 7(2/1) read as follows 

dyi g(S{x,t)) - (yi + y 2 ) dy 2 g(S(x,t))-y 2 , . 

" t e ' ~df = T a ' t = ( 3 - 2 ) 

where t e <C t a are constants, x is the current position of a cell, S : ffi™ x [0, 00) — > [0, 00) is the concentration of 
the chemoattractant, and g : [0, 00) — > [0, 00) models the first step of signal transduction. The constant 70 > is 
the turning rate if no chemoattractant is present, which is changed by the linear function —j3yi, with /3 > 0, if 
attractant gradients are present. 

In this paper, we restrict the random walks to movement along the real line, which means that individuals 
move to the left or to the right with constant speed s and, at random instants of time, change their direction 
with turning frequency 7. In this case, using l|3.2l) in suitable parameter regimes, one can derive a macroscopic 
partial differential equation for the density of individuals N = N(x, t) of the following form |1(J| . 

d 2 N dN d ( 2 8N 2J3s% \ 

The macroscopic equation (|3.3|l is valid for shallow gradients of the signal (small S'(x)) and for a suitable order 
of magnitude of the parameters involved (see ^3 for details). 

Since, bacteria are too small to sense spatial gradients of the chemoattractant over their body lengths, they 
alter their turning rates as described above, to achieve the desired response to changes in chemoattractant 
concentration. On the other hand, eukaryotic unicellular organisms like Dictyostelium discoideum are large enough 
to sense directly the chemical gradients and respond to them appropriately. Motivated by this observation, in 
the following section we present a simple example of a ID random walk of individuals such that a cell can sense 
directly the gradient of chemoattractant S'(x) and respond with changes of its direction according to the gradient 
seen by the cell. 



3.2 Chemotaxis with directional sensing 

We consider the random movement of individuals which reduce their probability of changing direction when 
moving in a favorable direction, e.g., in the direction of increasing attractant. We suppose as earlier that a 
particle moves along the x axis at a constant speed s, but that at random instants of time it reverses its direction 
according to a Poisson process with turning frequency 

1 = l0 ±bS'{x) (3.4) 

where b is a positive constant and the sign depends on the direction of the particle movement: plus for particles 
moving to the left, and minus for particles moving to the right. Let R(x, t) (resp.L(x, t)) be the density of particles 
at (x,t) which are moving to the right (resp. left): then R(x,t) and L{x,t) satisfy the equations 

^T + S 1T = -(7o - bS'(x))R + (70 + bS'(x))L, (3.5) 
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Fit Fit 

- s°-g = (70 - bS'(x))R - (70 + bS'(x))L. (3.6) 

Equations of this type have been studied by many authors, and for a discussion of previous work see |32l 121) . 

The density of particles at (x, t) is given by the sum N(x, t) = R(x, t)+L(x, t), and the flux is sR(x, t) — sL(x, t). 
We are primarily interested in the evolution of the macroscopic density N, and therefore we rewrite the equations 
1)3.5)1 and 1)3.6)1 as the equations for the variables N and J given by 

N = R + L, J = R-L o R= j" , L = — ^ — , (3.7) 

where J is a rescaled flux. Then adding and subtracting 1)3.5)1 and 1)3.6)1 . gives 

9N dJ n 

lH +S lTx= > (3 - 8) 
r) T BN 

^ + = ~ 2 ^foJ + 2bS'(x)N. (3.9) 

Thus the random walk can be described by the closed system of two equations l|3.8|l and l|3.9|l with given initial 
conditions iV(-,0) and J(-,0). 

Finally, assuming sufficient smoothness, we can convert 1)3.8)1 - (|3.9|l into a second order damped hyperbolic 
equation for N, namely 

d 2 N „ ON ^d 2 N „, d 



dt2 '^ = S ^-2 b s-(S>(x)N). (3.10) 

This is a hyperbolic version of the classical Keller-Segel equation |22II23| . Note that (|3.1U|) has the same structure 
as (|3.3|l . which can also be written as a system of two equations of the form (|3.8|l - (|3.9|l . Therefore, the system 
(|3.8|l - (|3.9|) can also be viewed as a macroscopic description of bacterial chemotaxis. 



3.3 Scaling and mathematical formulation of main problems 



If we consider the system (|3.8|l - (|3.9|) as a description of the collective movement of bacteria E. coli, then we 
can give biologically realistic values for the parameters s and 70. The speed of a bacterium is s ~ 10/im/sec and 
the turning frequency is 70 — 1 sec -1 . To nondimensionalize equations 1)3.8(1 - 1)3.9(1 . we choose the characteristic 
time scale To = Iq 1 ! we denote the characteristic space scale as Lo, and the characteristic concentration as No. 
Define 



sT 
Lo ' 



S'(x) = 



bS'(x)T 



N = 



N 



J = 



L ' No' 
Then the nondimensionalized equations 1)3.8)1 - 1)3.9)1 have the form 



J_ 

iVn 



t = 



t 



x 

To- 



on „dJ n 



% + s d Ji = -2J + 2S'(x)N, 
dt ox 



and to simplify notation, we drop the hats in 1)3.12)) and obtain the nondimensionalized system 

ON dJ 



(NJK 



dt 
dJ_ 
dt 



~ s ir = ° 

ox 

9N 
ox 



2S'(x)N 



(3.11) 
(3.12) 

(3.13) 
(3.14) 



Here we have one dimensionless parameter s and one dimensionless function S'(x), and we estimate the orders 
of them as follows. In a typical macroscopic bacterial experiment the characteristic length scale Lo is 1 cm, and 
since the characteristic time scale is Tq = j^ 1 — 1 sec, we have s ~ 10 -3 . If the characteristic length scale is 10 
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cm then s = 10 4 , and in either case s is a small parameter. A realistic choice of S'(x) must ensure that the 
turning rate (|3.4[1 is positive, i.e., \S'(x)\ < 1. Hence, we will assume throughout that 

s<l, \S'(x)\<l. (3.15) 

The system (NJ) is a linear hyperbolic system of two equations with nonconstant coefficients, which can be 
rewritten in diagonal form as a system of two equations for the right and left fluxes (cf. (|3.5|) and (|3.(i|) L Thus, 
(NJ) can be rewritten as 



(RL) 



-ft+ s fa=-l 1 - S '( X )} R + [1 + S '( X )] L ( 3 - 16 ) 
- s°-g = [1 - S'(x)]R - [1 + S'(x))L (3.17) 



We also know that the system (NJ) can be written as a single second order equation for N (compare with (|3.10|) L 
or as the following system for the variables N and U. 



(NU) 



In the following sections we study the system (NJ) or its equivalent formulations (RL) and (NU). We will restrict 
our computations to the finite interval [0,2] with no flux boundary conditions which, in the formulation (NJ), 
can be written in the form 

dN dN 
J(0,i) = — (0,t) = S'(0) = 0, and J(2,t) = — (2,t) = S'(2) = 0, for t > 0. (3.20) 

ox ox 

As indicated here we also impose no-flux boundary conditions on the signal. 

Finally, let us identify the dimensionless times of interest. The characteristic time scale was set as to the mean 
turning time, i.e., T = 1 sec, since that characterizes the microscopic dynamics, but the macroscopic times of 
interest in pattern formation experiments are several hours or days. From the mathematical point of view we are 
interested in the long term dynamics and steady states, and therefore we want to develop methods to compute 
the density profile N(x,t) for dimensionless times i > 1. 



3.4 Slow and fast variables and the slow manifold 



In this section we consider spatial regions where the signal derivative is either zero or maximal possible (to assure 
a nonnegative turning rate) . We show that in such regions the fluxes relax to functionals of the density for large 
times, i.e. the memory of the initial flux decays quickly. Thus the long-term dynamics can be described by a 
single first-order in time equation for the density N. Similar conclusions can be also made about systems (RL) 
and (NU). For example, in the case of (RL), we could characterize the long-term dynamics using only the right 
flux R. or only the left flux L, or any linear combination of R and L (e.g., the density N). Knowing the density 
AT, we can compute either (or both of) the right and left fluxes - alternatively, these fluxes quickly evolve to 
functionals of the density field; this constitutes our "slow manifold" . The choice of the "right" observables can 
be made by the modeler; for historical (as well as practical) reasons we will use the density N in the following as 
a description of the slow variables. 



3.4.1 Special choices of S'(x) 

If S"(x) = 0, then system (NJ) can be rewritten as a second order damped wave equation 

d 2 N ON 2 d 2 N 

-w +2 w= s ^- (3 - 21) 
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It is well-known 02| that the asymptotic behavior of the solution of (|3.2I|) under the boundary conditions (|3.2(JI) 
is given by the corresponding diffusion equation 

dN _ s 2 d 2 N 

dt - 2 dx 2 ■ {6 ] 

Consequently, the long-term, slow dynamics can be described by this first order in time equation for density only. 

Next consider a spatial region where the signal gradient is the maximum possible, i.e., S'{x) = 1. If the region 
with maximal signal gradient is large enough, then (RL) in this region reduces to 

dR OR dL dL , 

Seen the leftward flux L decays exponentially according to the second equation, the long-term behavior (in large 
spatial regions with S'(x) = I) is given by the rightward flux only. Since, N= R + L and L — > 0, the long time 
dynamics is simply described by the first order transport equation 

dN dN „ , nnJ . 

A similar transport equation holds for the minimal possible signal gradient S'(x) = — 1. Of course the boundary 
conditions (|3.20() require that we cannot choose S'(x) = 1 in the whole domain of interest, and consequently, 
(|3.24() only gives a good approximation of the behavior of cellular density in large spatial regions with maximal 
signal gradient. On the other hand, if we consider the random walk in a finite domain and we look for long term 
dynamics/stationary state then the no-flux boundary conditions H3.20J1 have to be taken into account. 



3.4.2 (NJ) for general signals 

For general signals, the behavior is just a combination of transport and diffusion as given by the second order 
equation H3.I0[1 . The steady state of (NJ) under no-flux boundary conditions is given by 

d 2 N d 

and it follows that 

N s (x) = Cexp (^S{x)^j (3.25) 

where the constant C is given by the initial condition for N . The interesting question is whether the behavior 
of (NJ) can indeed be described by a single first order equation for large times. The simplest choice is to use a 
parabolic counterpart of l|3.10fl . given in dimensionless form as 

dN s 2 d 2 N d 

Equation H3.26|l has the same steady state as (NJ), and moreover it reduces to (I3.22|l for constant signals. On the 
other hand, if S'(x) = 1, then equation 1)3. 26|) differs from Ij3.24|l by the term which adds artificical diffusion 



to the system Consequently, if we have extended spatial regions where S (x) — 1, then the equation 
gives different transient behavior than (NJ) but finally leads to the same steady state as (NJ). It is important to 
note that a major issue in equation-free computation is how many independent variables are needed in order to 
close with a first order in time system, because it may be difficult to initialize microscopic variables consistently 
with given macroscopic observables and their history {e.g. their first order time derivatives). In regimes where 
at least a second-order-in-time equation is needed for closure, initializing the density is not enough; the time 
derivative of density must also be prescribed. In such a case we would use an alternative initialization for 
equation-free computations: we would prescribe right- and left- going fluxes R and L, which would be sufficient 
to start a particle-based simulation, because it is much easier to initialize particles based on more than one 
independent variables rather than based on the "history" of a single variable. 2 

2 In doing projective integration based on simulations over the entire spatial domain, the spatial order of the equation does not 
play a crucial role. If, however, one tries to use equation-free techniques such as the gaptooth scheme and patch dynamics 24 11811351 . 
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Figure 4: (a) Graph of "hat-profile" signal function S(x). (b) Graph of S (x). (c) Graph of S (x). 



3.5 Test family of signal functions 

In later sections, several numerical computations are presented. Here we introduce the test family of signal 
functions which we will use in these illustrative examples. In all examples, we consider the problem (NJ) on 
the interval [0,2] with no-flux boundary conditions, where the signal belongs to the one-parameter test family of 
signal functions given by 

S a (x) = aS(x), for a e [0, 1], (3.27) 

where S(x) is a fixed signal function and a scales the strength of the signal. The signal function S(x) is chosen 
in the following form (see also Figure 0}: 
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Since the maximal absolute value of the derivative S (x) is equal to 1, the assumption l|3.15f) requires that 
a G [0,1], and a — 1 means that the signal derivative S' a (x) is maximal possible in some subintervals of the 
domain [0,2]. For the S(x) chosen, the signal gradient S' a {x) is zero in the intervals [0, g] , [|, |] and [§,2] , so 
the behavior will be similar to the diffusion equation there (for any a). If a = 1 in l|3.27fl , then the signal derivative 
S[(x) is maximal possible, equal to 1, in the interval [|, |] ; consequently, the right moving individuals will never 
turn in this interval and the corresponding coarse equation is a transport equation (|3.24|l there. Similarly, the 
signal gradient is minimal, equal to - 1, in the interval [|, |] ; consequently, the left moving individuals will never 
turn in this interval and the corresponding coarse equation is again the transport equation there. 



4 Projective integration 

The next objective is to study the so-called projective integration of the system (NJ), or its equivalent forms (RL) 
and (NU). To that end, we first summarize results from about the projective forward Euler method. Suppose 
that we want to solve the initial value problem for the linear system of ordinary differential equations 

^ = £y, 2/(0) =y 0) (4.1) 

where y is n— dimensional vector and C is a n x n matrix of real numbers. Given constants k and M and step 
size St, the projective forward Euler method (P^) can be described as follows |16j: 

implementing effective matching conditions between patches becomes important, and that is crucially affected by the spatial order of 
the effective evolution equation. The design of computational experiments to determine the spatial order of an unknown (in closed 
form) equation is an interesting subject, discussed in part in I2KI . 
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(Pjjr-1) Use the forward Euler method 3 to integrate the system (|4.1H over fc time steps 
of the length St to compute y(t + kSt) from y(t); 



(P k -2) perform one more integration step to compute y(t + kSt + St) from y(t + kSt); 

>A 



(P^ f -3) perform an extrapolation over M steps, using y(t + kSt + St) and y(t + kSt) to 



estimate y{t + kSt + St + MSt) as y{t + kSt + St + MSt) = (M + l)y(t + kSt + 
St) - My(t + kSt). 

Thus, the procedure (P^-l) - (P^-3) integrates the system over the (fc + 1 + M) steps of the length St. Next, 
we have the following result 

Lemma 1 Method (Pff -1) - (Pf? -3) for solving \4-l\ * s stable provided that the error amplification given by 

a(\St) = \(M + 1)(1 + XSt) - M~\ (1 + XSt) k (4.2) 
satisfies \a(X5t)\ < 1 for all X in the spectrum of the matrix C of system ]4-H - 



Proof: See ^1 where a more general linear stability analysis for systems of nonlinear ODEs is done. 

The absolute stability region in the complex A(5t-plane, which is plotted in Figure [S] and Figure Efa), is the area 
inside the curve |<r(A<5t)| = 1. We see that the region splits into two parts for large M. Consequently, the constant 
M can be large if the spectrum is concentrated into two widely-separated regions corresponding to the fast and 
slow components. We also see that we increase the part of the stability region corresponding to fast components 
if we increase the number of inner integration steps k. The stability region for fc = 1 is given in Figure |^a), for 
fc = 2 in Figure [SJa) and for fc = 10 in Figure GJb). 

In the following sections, we discretize the PDEs using the method of lines. Some of the systems we study will 
have a real-valued spectrum for parameter values of interest. Consequently, the interesting part of the stability 
region from Figure EJa) is its intersection with the real axis. For large M, the real stability region comprises the 
union of two intervals given by Lemma [5] for k — 1, and plotted in Figure EJb). 



(a) 



k=2 



M=5 
M=8 



CO . 

Q- i : 

>. ! 1 

a \ : 

.£ \ \\ 
05 \ \; 




-9 



M=10 
M=20 



(b) 



0.5 



-0.5 



-1 -0.5 
real part of A.dt 




-1 -0.5 
real part of Xdt 



Figure 5: (a) The regions of absolute stability of P k methods for k — 2 and M = 5 (dot-dashed line), M = 8 
(dotted line), M = 10 (dashed line) and M = 20 (solid line), (b) The regions of absolute stability of Fff methods 
for fc = 10 and M = 20 (dot-dashed line), M = 30 (dotted line), M = 40 (dashed line) and M = 60 (solid line). 



Lemma 2 Suppose that the eigenvalues of the matrix C are all real. Then the procedure (Pff-1) - (P^f -3) with 
fc = 1 and M > 5 for solving \4-l\ is stable provided that 

XSt G (C, B) n {A, 0), for all X in the spectrum of X, (4.3) 



3 In fact any other integration scheme can be used here. 
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Figure 6: (a) The regions of absolute stability of methods for k — 1 and M — 2 (dot-dashed line), M = 3 
(dotted line), M — 4 (dashed line) and M = 6 (solid line), (b) Intersection of stability region from part (a) with 
real axis plotted as as a function of M. The equations for boundary curves A(M), B{M) and C{M) are given in 
Lemma\^ 



where C < —1 < B < A < are given by 



C = -l-' B = _M + 2 + y/{M-2F-8 M + 2-V(M-2r-8 

M + 1 2(M + 1) ' 2(M + 1) v ; 



Proof: This is an easy consequence of LemmaQ] Q.E.D. 

From Figure Efb) we see that (in the case of real spectrum) one can choose a large projective jump M provided 
that the spectrum of C lies in two small intervals, separated by a spectral gap. Later, we will see such linear 
systems arising in our simulations; the natural question then is: if we know (or can estimate) the spectrum of C, 
what is the maximal possible choice of M such that the Pj^ method is stable? The answer is given in the following 
lemma. 

Lemma 3 Suppose that eigenvalues of the matrix C are all real. Let —2 <c<— l<6<a<0 and St > be 
given constants such that 



X5t 6 (c, b) n (a, 0), for all A in the spectrum of C. 



(4.5) 



Then (P k -1) - (P k -3) is stable for all M satisfying the inequality 



M < min 



l + (l + q) fc+1 
a(l +a) k 



l + (l + 6) fc+1 
b{l + b) k 



l-(l + c) fc+1 
c(l + c) k 



(4.6) 



Proof: The amplification factor (|4.2|) is given by the formula 

a(XSt) = MXSt(l + X5t) k + (1 + XSt) k+1 . 
In order to have a stable method, the following three inequalities must be satisfied simultaneously: 

a(a) > -1, a(b) > -1, and a(c) < 1. 

Solving for M, we obtain Q.E.D. 

Finally, let us note that the results of this section could be also viewed as results of linear stability analysis of 
projective integration of general nonlinear systems of ODEs of the form y' — F(y), y(0) — yo, where y is an 
n— dimensional vector and F : R" — > M™ [Ttj] , 
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4.1 Projective integration of chemotaxis systems (NJ), (RL) and (NU) 



Before we implement coarse projective integration, we illustrate the use of projective integration and the factors 
affecting its implementation and effectiveness through the use of discretizations of the chemotaxis equations 
themselves. In this context, it is convenient to think that we only have available as an "inner simulator" a black- 
box dynamic integrator with a small, fixed time step (for example, a forward Euler simulator of a discretization 
of the problem), and that we are attempting to accelerate this black box code. 

Given a signal profile S(x), the speed s and initial conditions, we will look for the solution of (NJ) in the 
finite interval [0,2] with no flux boundary conditions (|3.20[1 . To do that, we will discretize (NJ) and rewrite it 
as a system of ordinary differential equations of the form (|4.1() using the method of lines. The resulting system 
of ODEs is a starting point for our basic projective integration algorithm, which is little different than (Pj^-1) 
- (Pj^-3). It is based on the sketch in Figure El Choosing a suitable time step St, and constants k and M, 
the algorithm is given in the following three steps (Prf -1) - (Prf -3). Note that the steps (Prf -1) - (Pr^-2) 
correspond to the step (i) as outlined in Figure [21 and the last step (Pr^-3) corresponds to step (ii) in Figure [21 

(Pr^-l) integrate system (NJ) over k time steps of length St to compute N(t + kSt) and J(t + kSt) from 
N(t) and from suitably initialized flux J(t) - see jO} and (|4.10l) ; 

(Pr£ f -2) perform one more inner integration step to compute N(t + kSt + 5t) and J(t + kSt + St) from 
N(t + kSt) and J(t + k5t); 

(Pr^-3) perform an extrapolation over M steps, using N(t + kSt + St) and N(t + kSt) to compute N(t + 
kSt + St + MSt) = (M + l)N{t + kSt + St) - MN(t + kSt). 

Note that we can approximate the time derivative of N in step (Pr^-2) by 

dN _ N(t + kSt + St) - N(t + kSt) 
~dt ~ St 

and therefore the step (Pr^-3) is equivalent to 

dN 

N(t + kSt + St + MSt) = N{t + k5t + St) + MSt— 

which is the forward Euler projective step. We see that step (Pr£ f -3) is really equivalent to step (ii) in Figure 
121 It is important to notice that integrating the full system (NJ) requires initialization not only of the density N 
(which is prescribed) but also of the flux J, which is not; this will be discussed further below. As we mentioned 
in Section |21 the coarse/projective integration method is efficient provided that we can choose a large projective 
time T in step (d) in Figure [21 relative to the time At of the steps (a) - (c) from Figure |21 and still retain accuracy 
and stability. Using the notation from Section [21 we have 

At = kSt + St, T = MSt, (4.7) 

and consequently, the gain Q of the method (|2.2|) can be expressed as 

„ T + At M + k+1 



At k + 1 



(4.8) 



Our goal is to make this gain as large as possible. Moreover, in order to use the scheme (Prf -1) - (Prf -3), we 
have to specify the spatial discretization of (NJ). We study two options in Section |4~2*1 Finally, we also have to 
specify how we initialize the flux in step (Prjj^-1). There are several possibilities for doing this, the easiest of 
which is to use the initial flux J(t) in step (Pr^-l) given by 

J(t) = 0. (4.9) 

We can also use as an initial guess the value of the flux computed in the previous step (Pr^-2) corresponding to 
a time MSt ago, i.e., before the projective jump. Thus we could use 

J(t) = "flux J(t - MSt) which was computed in the previous step (Pr^-2)" (4.10) 

A more sophisticated flux initialization is used in Section 03 which deals with Monte Carlo simulations (see also 
Figure [3J. 
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4.2 Discretization of (RL) and (NU) 



Various possibilities exist for discretizing the system (NJ) in the spatial domain; we start with one which is 
based on the equivalent form (RL) and on upwinding. The advantage of upwinding is that it provides a more 
stable scheme for problems with a significant convection component, but it introduces artificial diffusion into the 
problem [30]. Another possibility to spatially discretize (NJ), (RL) or (NU) is to use central differences, which 
leads to equation (|4.1ti|) . 

First, to solve the system (NJ) numerically, we transform it to the system (RL) of two first order equations 
in diagonal form. We want to solve (RL) over the interval [0,2] with boundary conditions given by l|3.20[l . We 
choose a number n and a mesh size Sx = 2/n, and we discretize the interval [0,2] with n + 1 mesh points 



Xk — k ■ Sx, for k = 0, . 



, ,n. 



(4.11) 



Next, we define 



Ri{t) = R(xi,t), Li(t) = L(xi,t), and % = S'(xi), 



0. 



The zero flux boundary conditions l|3.20[l simply mean that Rq = Lq and R n = L n ; consequently, we have to 
compute the time evolution of the 2n— dimensional vector 



w 



(R%, i?2, • • • R n -i, Rn, Lo, L\, L2, ■ 
To discretize spatial derivatives in (RL), we use upwinding, that is, 



(4.12) 



dR 

dx 



(Xi,t) 



Ri(t) - flj-i(t) dL 
Sx dx 



(Xi,t) 



i(t)-Lj(t) 
Sx 



Then, the solution of (RL) with boundary conditions <|3.20|) is approximated by the solution of a system of 
ordinary differential equations 

—rr = Aw, w(U) = Wq, 



dt 

where wq is a given initial condition and matrix A is defined by 



A- 



e 





l-S[ 






-1-e + SS 






1-Sa 






e + S' n . 
e 






l + S[ 





e 

L-e- 




Si 





i + s' 2 





e 

-1-e-Sk 






l+S'n. 








(4.13) 
\ 



where 



s 

Sx 



(4.14) 



Since we have approximated the original PDE system as a system of ordinary differential equations of the form 
(|4.1(l . the results in Lemma Lemma and Lemma can be applied. Alternatively, we can discretize the 
chemotaxis system in its equivalent form (NU) using standard central differences to approximate spatial derivatives 
in (NU). We use n + 1 mesh points H4.11(l and we define 



Ni(t) = Nfat), Ui(t) = U(x it t), and S[ = S'{xi), 



0,.. 



, n, 



z = (N , Nx,N 2 , N n - 2 ,N n -!, N n , U , U u U 2 ,..., t/„_ 2 , £7 n _i, U n y 



(4.15) 
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where / is (n + 1) x (n + 1) identity matrix and e is given by (|4.14|) . Then, the solution of (NU) with boundary 
conditions i|3.2(Jfl is approximated by the solution of a system of ordinary differential equations 



dz 
dt 



Bz, 



z(0) = z , 



(4.16) 



where zq is a prescribed initial condition. 



4.3 Efficiency of projective integration 



First, suppose that there is no signal gradient in the domain of interest, i.e., we put S' Q = S[ = . . . S' n — in 
matrices A and B. Choosing n = 40, the real parts of the eigenvalues of A and B as a function of e are plotted 
in Figures[7|[a) andOJa), respectively. We see that there is a clear spectral gap for small e. The eigenvalues of A 
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Figure 7: (a) Graph of real parts of eigenvalues of matrix A for n — 40 and no signal gradient in the environment, 
i.e., S' = S[ — . . . S' n = 0. The eigenvalues are real for e € [0, 1] and there is a spectral gap between A41 and A42. 
(b) Graph of real parts of eigenvalues of matrix A for signal given by \3.27l with a = 1 and for n = 40. 



are all real for e £ [0, 1] and they satisfy 

A £ 

The eigenvalues of B are all real for e £ [0, 0.5] and for no signal in the environment and they satisfy 



-2e,0l |J (-2-25,-2 



A e 



-1 + ^1-462,0] (J [-2,-1- v 7 !-^ 2 



(4.17) 



(4.18) 



The spectral gap between —2e and —2 in the case of matrix A is independent of the signal as can be seen from 
Figure EJb) , where we use the signal profile (|3.27l) with a = 1. We see that some eigenvalues changed, but that 
the spectral gap between A41 and A42 survived. The imaginary parts of the eigenvalues do not grow significantly 
with a, and consequently the values of the real parts determine the stability of the scheme; we can use results 
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Figure 8: (a) Graph of real parts of eigenvalues of matrix B for n — 40 and no signal gradient in the environment, 
i.e., S' = S[ — . . . S' n = 0. The eigenvalues are real for e £ [0,0.5] and there is a spectral gap between A41 and 
A42. (b) Graph of real parts of eigenvalues of matrix B for signal given by \3.2'1\ with a = 1 and for n = 40. 



from Lemma0for matrix A and small e £ [0, 1). To do that, we specify the time step St. Since we want XSt close 
to -1 for eigenvalues corresponding to fast modes, we put 



St = 0.5 



(4.19) 



Considering our scaling (|3.11|l . we see that l|4.19[) means that St is equal to time 1/(270). Next, if k is at least 2, 
then the "component" of the stability region around —1 is more extended then its second component around 
(see Figure EJl. Consequently, using Lemma the size of interval containing the slow eigenvalues determines the 
gain Q of the method. Using l|4.17[) for the matrix A, we have 



Q = 



1 + (l-e) 



fe+1 



which is approximately 



for small e. 



(4.20) 



e(l-e) fe (fc + l)' rr ~~~' J e(k + l) 

Note that the gain Q, given by (|4.20|l . is independent of the signal strength a and it can be very large for small 
e. On the other hand, as we will see in Section [4.41 the choice of small e will decrease the accuracy of the upwind 
discretization A due to the strong artificial diffusion of the scheme. 

Next consider the matrix B. The real parts of its eigenvalues as functions of e are plotted in Figure |S] We 
see that the "boundary" eigenvalues —2, — 1 + \Jl — 4e 2 , — 1 + \/\ — 4e 2 , of B are signal independent. The 
eigenvalues are all real for e £ [0, 0.5] and for no signal in the environment. However, if we increase the signal 
strength a some eigenvalues become complex, as can be seen from Figure El where we plot the "slow" eigenvalues 
close to zero in the complex plane for n — 200, e = 0.1, and for different signal strengths. Choosing St by (|4.19|) 
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Figure 9: A plot of "slow" eigenvalues of matrix B for n = 200 and e = 0.1. We plotted only eigenvalues close to 
zero for different strength of the signal a from \3.2T\I , namely: (a) a = 0, (b) a = 0.5 and (c) a = 1. 

and k, we can (for small signals) apply the results of Lemma|31to compute maximal possible projective jump M 



1G 



and, hence, to compute the gain Q by 14. 8|) for small e and for small signals. Using (|4.18|) and Lemma for the 
matrix B for small signal gradients, we have 

1 + (1 - 4e 2 )(' £ + 1 )/ 2 1 

Q = , ■ , which is approximately for small e. (4.21) 

(-1 + VI - 4e 2 )(l - 4e 2 ) fc / 2 (fc + 1) y e 2 {k+l) y 1 

We see that discretization ((4.1611 results in a very large gain Q for small e and for small signal gradients. On 
the other hand, if we increase the signal gradients, then the result 1(4.21(1 is no longer true, because complex 
eigenvalues can appear outside the stability region (compare Figure [5] and Figure [5J. For example, we see from 
FigureHJc) that the slow eigenvalues lie in the complex interval [—0.02,0] x [O.li, O.li] for a = 1. Consequently, 
the absolute values of the imaginary parts of the eigenvalues are much larger than the absolute values of the real 
parts and the result 14.21(1 is not applicable for large signals. 



4.4 Accuracy of projective integration 

As we see in 1)4.20(1 . ((4.2111 and 1(4.14(1 . choosing a larger Sx will make e smaller and we will have a larger gain Q for 
the projective integration. On the other hand, a smaller Sx will increase the accuracy of the numerical method 
obtained by 1(4.13(1 or 1(4.16(1 . The right choice of Sx depends on the underlying signal. If we have signals with 
sharp second derivatives and if we want to capture the detailed transient behavior accurately, we have to use a 
sufficiently small Sx. However, if we want to make use of the spectral gaps 1(4.17(1 or ((4.18(1 . we must assure that 
s <C Sx to have e<l. 

Two types of errors arise in these computations: (1) the error between the projective integration of 1(4.13(1 or 
((4.16|) and the corresponding solutions of 1(4.13(1 or 1(4.16(1 . respectively; and (2) the error between solutions of 
((4.13|) or 1(4.16(1 and the exact solution of (NJ). The error in part (1) is sufficiently small as will be seen in Section 
14.51 it is easy to estimate this error here, since the exact solution of 1(4.13(1 . 1(4.16(1 or even of (NJ) can be found 
through careful, error-controlled computations. For microscopic simulations though, when the corresponding 
macroscopic equation is not known, estimating these errors becomes an important task; fortunately, numerical 
analysis techniques for on line a posteriori error estimates have been extensively developed for continuum prob- 
lems, and can be naturally incorporated in equation-free computation [ 121 . For example, comparing results of 
the same computation with half the projective time step can be used to estimate the error of the scheme and 
control projective time step selection; comparable techniques for adaptive spatial meshing can also be used. It 
is, however, important to note one "twist" to traditional a posteriori numerical error estimates: errors due to 
the estimation scheme, e.g., due to fluctuations in stochastic simulations; this can be controlled through variance 
reduction schemes, either by brute force computation of several replica simulations or possibly through biasing 
for variance reduction |29j . Beyond adaptive time steps, adaptive mesh sizes and possibly variance reduction, we 
will discuss at the end of the paper the adaptive check of the level at which a macroscopic description closes, i.e.. 
the number of macroscopic variables required, or the dimension of the "slow manifold" . 

We will now discuss errors of type (2), i.e., errors between solutions of ((4.13(1 or ((4.16(1 and the exact solution 
of (NJ). We can numerically estimate those errors by comparing the solution of 1(4.13(1 or ((4.16(1 for different Sx. 
Representative results can be found in Figure ITU1 where we used s — 0.0001, a signal given by 1(3. 27|) with a = 0.1 
and St given by 1(4.19(1 . As is well known, the upwinding 1(4.13(1 discretization introduces artificial diffusion to the 
problem which makes the solution more dependent on Sx. The central differences discretization ((4.1611 is more 
accurate here (40) . 

4.5 Numerical examples 

Here we present illustrative numerical results. In view of ((3.15(1 . we choose 

' Sx = 0.01, £ = T" = 0.01; (4.22) 



10000' Sx 

and we consider 201 mesh points ((4.11|1 in the interval [0,2]. The time step St is given by 1(4.19(1 and the initial 
condition is 

N(x,0) = l, J(x,0) = 0. (4.23) 
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Figure 10: (a) Graph of a solution of (RL) given by w = Aw at time t = 10 4 for different choices of Sx. We used 
s = 0.0001. signal given by j3.27[ with a = 0.1, St given by an d initial conditions \4-%3t - Consequently, 

Sx = 0.01, Sx = 0.001 and Sx = 0.0001 correspond to e = 0.01, e = 0.1, and e = 1, respectively, (b) Graph of a 
solution of (NU) given by z = Bz at time t = 10 4 for different choices of Sx. The parameters are the same as in 
(a). 



We know from Section EH"! that the discretization ((4.16(1 gives rise to a sufficiently accurate solution of (NJ) for 
the choice of parameters ((4.22(1 . so we start with the discretization 1(4.16(1 first. We learned in Figure El that we 
can have a large gain Q of using projective inegration of 1(4.16(1 if the signal gradient is small; consequently, we 
consider the signal ((3.27|) with a = 0.1. The numerical results for k = 1 and M = 398 are given in Figure ITT1 
Here the gain is Q = 200 using definition (|2.2() . In Figure ITTl we compare the solution of system (|4.16() with the 
projective integration of 1(4.16(1 . We see that the errors between the projective integration of the system of ordinary 
differential equations (|4. and the solution of 1(4.16(1 are small. Since the discretization 1(4.16(1 gives a reasonably 
accurate solution of (NJ), we can view also Figure [n] a s a plot of the exact solution of (NJ). Consequently, what 
we presented appears to be capable of significantly speeding up an explicit forward Euler method for small signal 
gradients (see also (161 1251 12*5] ). 

The second numerical example in this section is based on the upwind discretization 1(4.13(1 . We know from 
Section 14.41 that the discretization ((4.13(1 provides a less accurate solution of (NJ) than 1(4.16(1 for parameters 
(|4.22(l due to the artificial diffusion of the upwinding scheme. On the other hand, the gain Q of the projective 
integration method (|4.13|) is independent of the signal stregth a. Consequently, we will present here results for 
a = 1, i.e., when the signal is maximal possible. If we compare the results obtained by projective integration of 
(14.13(1 and the corresponding plots of solutions of l|4.13|l . we again obtain small errors (results not shown) similar 
to those in Figure 1111 This would again support the numerical results from |15j concerning the accuracy of 
projective integration of ordinary differential equations. Instead, we compare the results of projective integration 
for two different choices of e with an accurate solution of (NJ). We use either (|4.22|l or 

' 0.001, £ = ^-=0.1. (4.24) 



10000 ' ' Sx 

Moreover, we use the initial condition H4.23|l and St given by ((4.19(1 : the results are shown in Figure [El We see 
that the long time behavior is highly influenced by the artificial diffusion of the scheme. Projective integration 
with small e has large gain Q, but it will reach the steady state much faster than the exact solution of (NJ). Note, 
that it is not an inaccuracy in projective integration per se; the inaccuracy is created by the inaccurate spatial 
discretization, based on upwinding with small e. 

Next, we will discuss how the ideas described so far in this paper can be used in Monte Carlo simulations of 
chemotaxis. 



18 



t=0 



in 

ro 1.5 

"O 

> 

~ 1- 

o 
>. 

m 

5 0.5 

T3 



0.5 



1 

position 



1.5 



1-6 t=10 4 





1.06r 


W5 


1 .04 


CO 




Z3 






1.02 


> 










1 ■ 




"o 








c/5 


0.98 


£Z 




(D 




"O 


0.96 











t=1000 



0.5 



1 

position 



1.5 




t=10 5 















0.5 1 

position 



1.5 



t=10' 



0.5 




1.5 



position 



position 



Figure 11: The time evolution of the density of individuals for s = 0.0001 and a = 0.1. We plot the solution of 
the system \4-lb\ (dashed line) and the solution obtained by the projective algorithm (P fe -1) - (Pff-3) for \4-lb\ 
with k = l and M = 398 (solid line). The gain Q is 200. We use (fcZQ , (gT^ and $Jfy . 
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Figure 12: (a) Density of individuals for s = 0.0001 and a = 1 at time t = 1000. We plot the solution of (NJ) 
given by accurate numerical method (solid line), the solution obtained by projective algorithm (P^ f -1) - (Pj^-3) 
for \4- with \4-2'£\ , k = 1 and M = 198, i.e., with Q = 100 (dashed line). We also plot the solution obtained 
by projective algorithm (P^-l) - (Pff-3) for \4-l?l[ with \4-24\ , k = 1 and M = 18, i.e., with Q = 10 (dot-dashed 
line). In all computations, we use \4- l&t an d \4-%ty - (b) The same plots as in (a) for time t = 10 5 . 



5 Coarse projective integration: a kinetic Monte Carlo example 

In the previous sections we studied the gain Q of projective integration for the system (NJ) of deterministic partial 
differential equations. Here we present results of Monte Carlo simulations of the underlying random walks using 
coarse projective integration |15U17ll2"4"| that will make use of the previous analysis. While in principle we simulate 
the evolution of the particle density profile over the entire spatial domain, we will demonstrate how to perform 
the computations required for coarse projective integration on a relatively small portion of the full domain. This 
is based on the presumed smoothness in physical space of the evolving density profile, which constitutes the 
underpinning of equation- free methods such as the gap-tooth scheme |18l 1351 124| as described below. Here we are 
able to speed up the kinetic Monte Carlo simulation about a thousand times. 

Suppose that we have 2no random walkers in the interval [0,2], and suppose that we have only a kinetic 
Monte Carlo simulator to model the evolution of the system. As before, the interesting macroscopic quantity 
is the density of random walkers N which can be obtained as follows. We choose a macroscopic mesh size 6x 
and we discretize the interval [0,2] using mesh l|4.11|l . Then we obtain the (probability) density N i+ i/ 2 (t) at 
point ' +1 as the number of particles in the interval [&,, a^+i] divided by uqSx. We thus create a histogram of 
particles, which can of course be noisy. 

If we have randomly walking noninteracting particles, the histograms obtained by using 10 6 or 10 9 random 
walkers appear roughly the same; the former is just less "noisy" than the latter. Consequently, we can obtain 
relatively accurate results quickly by simply decreasing the number of particles. However, in many interesting 
biological problems, cells change their environment, they consume nutrients, secrete waste, etc. Consequently, 
cells interact through environmental chemicals and then the number of cells is prescribed by the biological setup, 
and we cannot change it without changing the computed solution. 

Therefore, in the examples of this section we will suppose that we do not know that the particles are noninter- 
acting; we will suppose that there is a fixed number 2uq of individuals in the domain of interest - the interval [0,2] 
- which are moving according to the rules of the random described in Section 13.21 We will show that in the case 
of a fixed number of particles the coarse integration method leads to an even larger gain Q than the projective 
integration method used earlier, where as before, the gain Q is defined by H2.1(l . In the numerical example, we 
choose 

2n = 10 8 , * = ^5> fo = 0.01, (5.1) 
i.e., we consider 201 equi-spaced mesh points (14. 1 If) in the interval [0,2] on which we track the evolution of the 
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macroscopic density, and the parameters are the same as in (I4.22|) . 

Monte Carlo simulations are performed as follows. Each particle is described by two variables - position 
x € [0,2] and velocity ±s. We use a small microscopic time step dt — 0.01, i.e., the unbiased turning frequency 
divided by 100, and during each time step the particle moves with speed s in the chosen direction. At the end of 
each time step, a random number chosen from a uniform distribution on [0, 1] is generated and compared with the 
probability of the turn jdt = (1 ± S'(x))dt. If a turn occurs, the cell will move in the opposite direction during 
the next time step. To apply the previous results, we choose a macroscopic time step St given by (|4.19|) and do 
kinetic Monte Carlo simulations in the interval [t, t + St], which means that we use the Monte Carlo simulator for 
5t/dt microscopic time steps dt. 

Since the histograms are noisy, we will work with the integral of the density - i.e., with the cumulative density 
function defined by 



C{x,t) = / N(x,t)dx. (5.2) 
Jo 

Discretizing the interval [0,2] using mesh (|4. 1 1|> - we obtain 

i 

C i {t)=C(x i ,t) = 5x^2N k _ 1/2 (t), and C„(t) - C^t) = ^_ 1/2 (t)fa. (5.3) 

fc=i 

In particular, the number of particles in the interval [0, Xi], i = 1, . . . , n is given by 

i 

n Ci{t) = n Q 8x^N k _ l/2 (t), i = l,..,,n. (5.4) 

k=l 

In order to use coarse integration, it is important to compute the change of Ci(t) during the time interval [t, t + 5t]\ 
equivalently, we want to know the change of the number of particles in [0, Xi] during the time step [t , t + St] . Given 
that the speed of the particles is s, only particles which are in the small interval [xi — sSt, Xi + sdt] at time t can 
enter or leave the interval [0, xi\. Consequently, only a small number of particles around each mesh point have 
to be simulated (compare Figure IT^T b)): of course we are implicitly assuming that the discretization mesh is fine 
enough so that interpolation between mesh points provides an accurate estimate of the evolving density profile. 
Using H4.19fl and previous results, we choose 

<5t = 0.5, T = 99, (5.5) 

and compute the cumulative density at time t + 2St + T = t + 100 from the cumulative density function at time 
t by the following algorithm (compare with Figure UHIa) and Figure 0) 

(al) Given a macroscopic initial cumulative density C(t) at mesh points l|4.11|l . we compute the density 
Ni-i/2(t) by the formula 

d(t) - d-^t) . 

Ni-i/ 2 {t) = ^ , i = l,...,n. 

We put noN i _ 1 / 2 {t)5x particles in each interval [xi-\,Xi] and distribute them so that the resulting prob- 



ability density function is a continuous piecewise linear function with value A r i _ 1 / 2 (i) at point 
i= 1, . . . , n. Thus (see Figure [r^T bl'l 



Xi — \-\-Xi 



N(x,t) — N i _ 1 / 2 (t) H ' — ' \x ) forxe 



•Ei — 1 Xi -f" Xi-^-i 



Moreover, we assign alternating velocities to the particles, so that the initial flux is effectively zero. As we 
mentioned earlier, we do not have to simulate all particles in [ z '~ 2 * , z f'" 1 " 1 ]; instead, we consider only 
particles which are inside a small interval [xj — 2sSt, Xi + 2sSt] around the macroscopic mesh point x^ (this 
could be thought as analogous to the gap-tooth scheme JH!> except that one does not have to formulate 
and impose effective smoothness boundary conditions, see Figure IT^T b)). 
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Figure 13: (a) Schematic of steps (al) - (d) of the coarse integration algorithm. Monte Carlo simulation is 
denoted by solid lines. Dashed lines denote relatively fast steps, i.e., resetting the values of density to N(t) in 
steps (a2) and (a3) and the projective step (d). The dot-dashed line represents the evolution on the slow manifold 
of the system. We assume that the slow manifold can be accurately parametrized by density, (b) Schematic of 
three macroscopic mesh points Xi-\, Xi and Xi+\ and kMC computational domains around them. Only particles 
close to the mesh points need be considered; the remaining particles will not leave/enter the interval [0,Xi] during 
steps (al) - (b2) and consequently, they do not have to be simulated. At the top, is the (piecewise linear) estimated 
density profile which is used in step (al). We place particles in the small computational domains that their number 
is consistent with this density profile. 



(a2) Evolve the system using the microscopic Monte Carlo simulator for time St. Then return the particles 
to their initial position as given in (al) but with a velocity equal to the values computed in (a2) (this and 
the following are preparatory steps to bring the microscopic initialization close to the slow manifold). 

(a3) Repeat (a2) again, i.e., evolve the system using the microscopic Monte Carlo simulator for time St. 
Then return the position of particles to their initial values as given in (al) keeping the velocities equal to 
computed velocities in (a3) (this can be repeated a few times). 

(bl) Using the positions and velocities produced at the end of step (a3), evolve the system using the 
microscopic Monte Carlo simulator for time St. Compute the number of particles in the interval [0, Xi] at 
time t + St for i — 1,2, ... ,n. 

(b2) Evolve the system using the microscopic Monte Carlo simulator for another time step St. Compute 
the number of particles in the interval [0, Xi] at time t + 2St for i = 1, 2, . . . , n. 

(c) Using data from (bl) and (b2), compute cumulative densities C(t + St) and C(t + 2St) at mesh points 
xq, xx, . . . , x n (this is the restriction step in equation- free computation). 

(d) Estimate the time derivative 

3d Ci{t + 2St)-d{t + St) 



dt ' St 
and take an extrapolation (projective) step 



(5.6) 



C l {t + 25t + T) = C l (t + 25t)+T^. (5.7) 

Then use Ct(t + 2St + T) as the new initial condition in step (al). 

The steps (al)-(d) of the algorithm are illustrated in Figure IT3T a^ where the slow manifold in density-flux space 
is shown as a dot-dashed line. Note that the steps (al) - (a3) correspond to the step (a) from Figure[3 They are 

4 This step annihilates the correlations between the present and the initial velocities of a particle, and at the macroscopic level the 
time required is essentially that in which a hyperbolic equation rather a parabolic equation is needed at the macroscopic level. This 
was already known to Einstein (c/. |.'S3I I. 
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preparatory steps used to initialize the flux close to the slow manifold (since we assume that the flux equilibrates 
quickly); they qualitatively correspond to evolving the macroscopic PDE for a short time constraining the density 
profile to be the one we want to prescribe as our macroscopic initial condition. Such constrained evolution 
preparatory procedures (like "umbrella sampling" ) are standard in computational chemistry [411 I34| . A more 
detailed description of such initialization algorithms in the case of legacy simulators can be found in (201 Q3| ■ 

The steps (bl) - (b2) correspond to step (b) from Figure Moreover, (bl) corresponds to the step (Prj^-1) 
and (b2) to the step (Pr fc -2) from projective integration algorithm of (NJ). Similarly, steps (c) and (d) can be also 
found in Figure[2 moreover, steps (c) and (d) together form step (Pr^-3) of the projective integration algorithm 
of (NJ). 

If there is a small number of cells in one of the computational domains, then the straightforward application 
of the algorithm (al) - (d) could give unrealistic results. For example, suppose that there are only two cells in the 

interval [aii i , atf] at time i, that the first cell moves to the interval [0, aJ»— i] during the time interval [t,t + 25t], 

the second cell moves outside the interval [0, Xj] and that no other cell crosses mesh points Xi—i and x% during 
time interval [t, t + 2St] . Then the time derivative of the cumulative density function l|5.6|) would be negative at 
point Xi and positive at Xi-\. Moreover, the projected solution (|5.7(l satisfies C;_i(i + 25t + T) > Ci(t + 25t + T); 
consequently there is a negative number of particles in the interval [aJj-i, Xi] at time t + 2St + T. To avoid this 
problem we have to consider more realizations for each computational domain containing a small number of 
particles, and compute an average over this ensemble of realizations. Practically, if the number of particles n, 
in the small computational domain around Xi is less than a given number m, we choose to repeat (al) - (d) for 
m/tii microscopic realizations in this computational subdomain. 

Numerical results for St = 0.5, T = 99, signal strength a = 0.1 and m = 10000 are given in Figure H~4l There 
are two sources of gain for this method. First, we have the gain of the projective step. In one step (al) - (d), 
we compute the evolution of the system over time T + 2St = 100 and we run the Monte Carlo simulator for 
time iSt = 2 in steps (a2) - (b2). Consequently, the gain factor of the projective step is (T + 2St)/4St = 50. 
The second part of the gain comes from the fact that important particles (for the estimation of the evolution 
of a smooth macroscopic density) are only those particles which are leaving/entering the interval [0,1,] at the 
endpoint. From Figure IT3T b). we see that only particles which are at time t with distance less than 2s8t = 0.0001 
from the endpoint can leave/enter the interval [0, x<] during steps (bl) - (b2). Consequently, only the fraction 
2no ^j^L — 4^2- of particles have to be simulated, and another factor of 50 appears in the gain. 

Therefore, the combined gain of the coarse integration and reduced spatial simulation (based on macroscopic 
density smoothness) is 50 x 50 = 2500. However, 2500 is not the actual gain Q because some computational time 
was lost by considering multiple microscopic realizations of domains which contained a small number of particles. 
In any case, we add less than 199771 particles to the simulation where 199 is the number of "inner" computational 
domains and m — 10000 is the minimal number of particles in each of them. Consequently, we actually simulated 
more cells than = 2 ■ 10 6 but, at any time, the number of simulated cells did not exceed + 199Af ~ 4 • 10 6 . 
So, in the worst possible case, we slow down the computation by a factor of 2, which means that the total gain 
G of the method is at least Q = 50 x 50/2 = 1250. 

In Figure [T4l we present the time evolution of the solution given by method (al) - (d) (solid line) compared to 
the solution of the macroscopic PDE equations (dashed line). Since the algorithm (al)-(d) computes cumulative 
density functions C{t) and we visualize the density N(t) in Figure ITU the results are noisy and the plots depend 
on the formula which is used to generate the density curves from the computed cumulative density data. To be 
precise, in Figure Il4l we show a plot of the function 

N( Xh t) = Ci+li \- x Ci - 1 ® . (5.8) 
Another representation of the results is given in Figure IT5l where we show results for time t = 10 4 using different 
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Figure 14: TTie time evolution of the density of individuals for (15. 1\ . (|5.5|l . a = 0.1 and initial conditions \4-2S[ . 
We plot the density given by coarse integration (al) - (d) and obtained by formula (|5.ff[l from the computed 
cumulative density function C{t) (solid line). Here we have gain Q = 1250. We also plot the solution of the 
corresponding macroscopic moment equations (dashed line). 
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(b) 




(c) 




position 



Figure 15: Plots of the density of individuals at time t = 10 4 for the same cumulative distribution function as 
in Figure \14\ using different expressions for for computing the discretized density, (a) formula (j5.9(l ; (b) formula 
(|5.1U|) ; (c) formula (|5.11|) . We also plot the solution of the corresponding macroscopic moment equations ( dashed 
line). 



formulas for the density function TV, namely 



(a) N( Xll t) = 



C i+2 (t) - Cj- 2 (t) 
A8x 



(b) N ( Xi+Xi + 1 t \ - C '«'+i(*) _ 



(c) N( Xi ,t) 



2 '7 Sx 

C l+2 (t) + C t+1 (t) - d-ijt) - Cj- 2 (t) 
6Sx 



(5.9) 
(5.10) 

(5.11) 



Comparing plots in Figure 1151 and the corresponding plot from Figure 1141 we see how the visualization of the 
results depends on the formula for estimating the discretized density function N(t) from the discretized cumulative 
density function C(t). In particular, if we use (|5.11(l instead of (|5.8(l in Figure ITfl then the results will look 
less noisy. Alternative discretizations of the particle density (namely, the use of orthogonal polynomials to 
represent the inverse cumulative distribution function, or ICDF) can be found in the literature ^JESj. Techniques 
for estimating smooth field profiles from noisy particle data have been proposed, among other places, in the 
computational materials science literature {e.g., the thermodynamic field estimator |27p. 

One can also decrease the noise in the computation, and the resulting macroscopic field estimates, by consid- 
ering multiple microscopic realizations, i.e., by increasing the value of m. However, if we increase m, then the 
gain Q will decrease (obviously the "wall clock" time of the overall computation remains the same if one does 
these computations in parallel). 

Certainly, there is a relationship between the initial number of particles 2no, the minimal number of particles 
in each small computational domain m and the gain Q of the method. If we have a large number of particles 2no, 
then we can choose a large m without significantly decreasing the gain of the method. On the other hand, if we 
have a stochastic problem with a small number of particles 2n , then it may not be appropriate to consider a 
closed PDE as a good model of a single system realization. In our example, m was chosen in such a way that the 
gain Q of the method was decreased only by a factor of two, and thus the Monte Carlo simulation was accelerated 
by a factor of at least 1250. Increasing m would further decrease the gain Q and reduce the magnitude of the 
fluctuations. 



6 Discussion 



In Section[5]we analyzed an example in which a simple coarse integration scheme was "wrapped around" a kinetic 
Monte Carlo simulation. The short (in time) bursts of kMC simulation were performed over only part of the full 
computational domain; this provides another important factor in decreasing the overall computational cost for 
such complex problems. The idea of reduced spatial as well as temporal simulation (the so called "gap-tooth" 
scheme and its combination with projective integration in "patch dynamics") is based on smoothness in the 
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evolution of macroscopic observables and constitutes a hallmark of equation-free computation. Let us note that 
the computation of long term dynamics of our system took several days on a IBM SP 375MHz Power3 processor 
using algorithm (al) - (d). Consequently, a computation using the kinetic Monte Carlo simulator would take 
several years and was not attempted. We estimated the accuracy of coarse projective computations by comparing 
to solutions of accurate macroscopic partial differential equations, which in this example happened to be known. 
When we do not have population level equations, we must use standard a posteriori error estimates to check 
accuracy and adaptively control the error of our results as discussed below. 

As we saw in Figure [3] for matrix B, the length of the possible projective step T, as determined by stability 
considerations, decreases with increasing strength of the signal a. The same is true for algorithm (al) - (d). If 
we increase a, then wc have to decrease T in order to have a stable scheme. In order to achieve stability for larger 
T we could use a similar strategy to that used for the matrix A: we could introduce artificial diffusion into the 
scheme which would make the scheme stable, independently of the strength of the signal a. It is not difficult to 
design a coarse integration scheme with artificial diffusion present; however, such a scheme would predict incorrect 
dynamics for the system. 

A better solution to the problem of large signal gradients is to note that large signal gradients are typically 
localised only in small parts of the domain of interest. In fact, the problem with the coarse integration scheme 
begins when a large signal gradient is present and particles become highly localized in space. Then the mesh is 
not fine enough in certain small domains of interest (around peaks) but it is sufficiently fine in the remainder of 
the interval [0, 2]. Similarly, the projective step is good for most of computational subdomains, but it would lead 
to instabilities because of strong signal gradients for a few of the computational domains. One could conceivably 
adapt the mesh, leading to a nonuniform mesh, finer in regimes with large signal gradients and coarser otherwise. 
Then we may need to make different projective jumps in different parts of the domain of interest; issues of this 
nature have been studied for nonuniform meshes in traditional continuum numerical analysis using adaptive mesh 
refinement (AMR) methods 0] , and in hybrid situations AMAR methods ^Sj • Efficient implementations of such 
adaptive techniques may be the key to significant acceleration of our illustrative Monte Carlo scheme, since they 
would allow us to obtain relatively accurate results for even larger sets of signal functions and for problems where 
the signal is also altered by the cells. 

Detailed methods have been developed for adapting the computation to the time and space scales of the 
problem in continuum numerical analysis. Adaptive stepsize selection in numerical integration, as well as adaptive 
mesh refinement in spatial discretizations is an indispensible part of modern software, and is typically based on a 
posteriori error estimates of the solution accuracy computed on line. These methods can be naturally incorporated 
in equation-free algorithms to control, for example, projective integration time steps to control accuracy. It 
should be noted that in addition to adaptive time-step selection (for coarse projective integration) and adaptive 
mesh selection (for gap-tooth algorithms), there is an additional type of adaptivity that arises in equation- free 
computation. This is the adaptive detection of the level of modeling, which may involve augmenting or decreasing 
the number of variables needed for closure. At a very qualitative level, adaptation of this "level of description" 
comes from the estimation of the gap between "fast" and "slow" system variables, which can be attempted using 
matrix-free eigensolvers. By initializing the microscopic distribution using more variables than the current level 
of modeling, one can try to estimate the characteristic relaxation times of the additional variables to functionals 
of the ones we need. This allows one to detect (while the level of description is still successful) whether variables 
that are treated as "fast" are becoming "slow" , and should be included as independent variables in the modeling. 
A good illustration of this is the evolution of stresses in a microscopic simulator of fluid flow: for a Newtonian 
fluid stresses rapidly become proportional to velocity gradients, while in non-Newtonian fluids this is not true, 
and one must use more independent variables to model such flows. This could be considered analogous to closing 
bacterial chemotaxis equations with only a single field (density) which can be done for long time dynamics in 
some parameter regimes versus needing two independent variables (right- and left- fluxes) to successfully close 
system in some other cases. In our case, the flux quickly becomes functional of density, as can be seen directly 
from simulations. 

A summary of the steps of our computational approach is as follows. 

• identify the appropriate level of closure 

• apply the equation-free computational algorithm 

• do a posteriori error estimation 
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As we discussed above, we have to first identify the level of closure, i.e. identify the slow dynamics of the 
system which we want to model. Then we can do coarse projective integration by making use of the spectral 
gap between fast and slow modes of the system. As we saw, it can be natural or desirable to combine coarse 
projective integration with gap-tooth methods, i.e. exploit the smoothness in physical space to only perform the 
computations on relatively small subdomains. As a result one can do transient calculations much faster than by 
direct simulations. If a modeller is interested in steady states and the transient dynamics are unimportant, then 
he or she can use other computational equation-free techniques (such as application of Newton-GMRES method) 
to obtain steady state behavior faster or do even bifurcation analysis |17l I87| . The final step is a posteriori error 
analysis as suggested above. This is an important issue if one wants to use our computational approach for the 
problems where macroscopic equations are unavailable. 

As we discussed, the large gain of the coarse projective integration is governed by the large spectral gap 
between fast and slow eigenvalues of the system. Our biological model system had such a large spectral gap 
because the mean running distance of individuals was much smaller than the size of the domain of interest. The 
method has a potential to speed up other models of biological dispersal with similar properties. 
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